/* 2010: Copyright by Massimo Bernaschi, Giorgio Parisi, Leonardo Parisi */
#include "cudamacro.h"
#include "mpiexch.h"

#include "hsgAMGPUGM.h"
#include "timing.h"

#include "dictionary.h"
#include "iniparser.h"
#include "util.h"
#if defined(DEBUG)
#include "util/cuPrintf.cu"
#endif

#include "MersenneTwister.h"
#if !defined(HOSTREDU)
#include "threadFenceReduction_kernel.cu"
#endif

extern void StartMpi(int *, int *, int *, char ***);
extern struct mpiexch *SetUpMpi(int, int, int);
extern void ExchMpi(struct mpiexch*);
extern REALENE SumEnergy(REALENE, int);
extern void StopMpi(void);
//-------------------------------------------------------------------------------------------------------
// Usage
//-------------------------------------------------------------------------------------------------------
void Usage(char *program) {
	fprintf(stderr,"Usage: %s -i inputfile [-r spin_file] [-o (over-relax only)] [-b (heat bath only)]\n",program);
}


//-------------------------------------------------------------------------------------------------------
// Init
//-------------------------------------------------------------------------------------------------------
void Init(void){
	
	if(sizeof(REAL)==sizeof(float)) {
		zero=0.0f;
		one=1.0f;
		two=2.0f;
		three=3.0f;
		zerofive=0.5f;
		zerotwentyfive=0.25f;
	} else if(sizeof(REAL)==sizeof(double)) {
		zero=0.0;
		one=1.0;
		two=2.0;
		three=3.0;
		zerofive=0.5;
		zerotwentyfive=0.25;
	} else {
		writelog(TRUE,APPLICATION_RC,"Invalid size of REAL: %d",sizeof(REAL));
	}
	
	zeroEne = 0.0;
	
	h_Ra = (REAL*) Malloc(ms_s);
	h_Rb = (REAL*) Malloc(ms_s);
	h_Rc = (REAL*) Malloc(ms_s);
	h_Ba = (REAL*) Malloc(ms_s);
	h_Bb = (REAL*) Malloc(ms_s);
	h_Bc = (REAL*) Malloc(ms_s);
	
	h_Jpx = (REAL*) Malloc(ms_j);
	h_Jpy = (REAL*) Malloc(ms_j);
	h_Jpz = (REAL*) Malloc(ms_j);
	h_Jmx = (REAL*) Malloc(ms_j);
	h_Jmy = (REAL*) Malloc(ms_j);
	h_Jmz = (REAL*) Malloc(ms_j);
	
	h_Energy = (REAL*) Malloc(ms_s);
	
	h_sdbL = findlog2(h_side);
	h_sdbLs = findlog2(h_side*h_side);
	
	h_sdbLm1 = h_sdbL-1;  //log_2(L) -1 used to divide by L/2
	h_sdbLsm1 = h_sdbLs-1;  //log_2(L*L) -1 used to divide by (L*L)/2
	h_maskL = (h_side>>1)-1;
	h_maskL2 = ((h_side*h_side)>>1)-1;
	h_maskL3 = (V(h_side)>>1)-1;
	h_Ldb2 = h_side>>1;
	h_Lsdb2 = (h_side*h_side)>>1;
	
	h_nspin=h_side*h_side*z_side/2;
	
	srand48(seed);
	
}


//-------------------------------------------------------------------------------------------------------
// CopyConstOnDevice
//-------------------------------------------------------------------------------------------------------
void CopyConstOnDevice(void){
	
	int tmp[2];
	
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_side",&h_side,sizeof(int),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_nspin",&h_nspin,sizeof(int),0,cudaMemcpyHostToDevice) );
	
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_sdbL",&h_sdbL,sizeof(int),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_sdbLs",&h_sdbLs,sizeof(int),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_sdbLm1",&h_sdbLm1,sizeof(int),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_sdbLsm1",&h_sdbLsm1,sizeof(int),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_maskL",&h_maskL,sizeof(int),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_maskL2",&h_maskL2,sizeof(int),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_maskL3",&h_maskL3,sizeof(int),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_Ldb2",&h_Ldb2,sizeof(int),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_Lsdb2",&h_Lsdb2,sizeof(int),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_h",&h_h,sizeof(REAL),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_Beta",&h_Beta,sizeof(REAL),0,cudaMemcpyHostToDevice) );
	
	tmp[0]=0;
	tmp[1]=-1;
	
	MY_CUDA_CHECK( cudaMemcpyToSymbol("jindred",tmp,sizeof(tmp),0,cudaMemcpyHostToDevice) );
	
	tmp[0]=-1;
	tmp[1]=0;
	
	MY_CUDA_CHECK( cudaMemcpyToSymbol("jindblue",tmp,sizeof(tmp),0,cudaMemcpyHostToDevice) );
	
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_const_one",&one,sizeof(REAL),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_const_two",&two,sizeof(REAL),0,cudaMemcpyHostToDevice) );
	MY_CUDA_CHECK( cudaMemcpyToSymbol("d_const_zerotwentyfive",&zerotwentyfive,sizeof(REAL),0,cudaMemcpyHostToDevice) );
	
	
	
}


//-------------------------------------------------------------------------------------------------------
// findlog2
//-------------------------------------------------------------------------------------------------------
int findlog2(int n){
	
	int i=0;
	
	if(n<1 || n&1) { writelog(TRUE,FINDLOG_RC,"Invalid value %d. Must be > 0 and even\n",n); }
	
	while(n>0) {
		n>>=1;
		if(n>0) i++;
	}
	
	return i;
}


//-------------------------------------------------------------------------------------------------------
// nextPow2
//-------------------------------------------------------------------------------------------------------
unsigned int nextPow2( unsigned int x ) {
	--x;
	x |= x >> 1;
	x |= x >> 2;
	x |= x >> 4;
	x |= x >> 8;
	x |= x >> 16;
	return ++x;
}


//-------------------------------------------------------------------------------------------------------
// FreeDevice
//-------------------------------------------------------------------------------------------------------
void FreeDevice(void){
	
	cudaFree(d_Ra);
	cudaFree(d_Rb);
	cudaFree(d_Rc);
	cudaFree(d_Ba);
	cudaFree(d_Bb);
	cudaFree(d_Bc);
	
	cudaFree(d_Jpx);
	cudaFree(d_Jpy);
	cudaFree(d_Jpz);
	cudaFree(d_Jmx);
	cudaFree(d_Jmy);
	cudaFree(d_Jmz);
	
	cudaFree(d_Energy);
	
}


//-------------------------------------------------------------------------------------------------------
// FreeHost
//-------------------------------------------------------------------------------------------------------
void FreeHost(void){
	
	Free(h_Ra);
	Free(h_Rb);
	Free(h_Rc);
	Free(h_Ba);
	Free(h_Bb);
	Free(h_Bc);
	
	Free(h_Jpx);
	Free(h_Jpy);
	Free(h_Jpz);
	Free(h_Jmx);
	Free(h_Jmy);
	Free(h_Jmz);
	
	Free(h_Energy);
	
}


//-------------------------------------------------------------------------------------------------------
// Dump_Spin
//-------------------------------------------------------------------------------------------------------
void Dump_Spin(int count){
	
	FILE* tempfile;
	char namefile[MAXSTRLEN];
	
	return;
	memset(namefile,'\0',MAXSTRLEN); // NON CONTROLLO SE TORNA BENE O MENO
	
	if(count>=0){
		
		snprintf(namefile,sizeof(namefile),"spin_dump_%d_%d",count,mpiid);
		
		MY_CUDA_CHECK( cudaMemcpy( h_Ra, d_Ra+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
		MY_CUDA_CHECK( cudaMemcpy( h_Rb, d_Rb+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
		MY_CUDA_CHECK( cudaMemcpy( h_Rc, d_Rc+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
		
		MY_CUDA_CHECK( cudaMemcpy( h_Ba, d_Ba+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
		MY_CUDA_CHECK( cudaMemcpy( h_Bb, d_Bb+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
		MY_CUDA_CHECK( cudaMemcpy( h_Bc, d_Bc+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
		
	} else {
		snprintf(namefile,sizeof(namefile),"spin_dump_start_%d",mpiid);
	}
	
	tempfile = Fopen(namefile,"w");
	
	for(int i = 0; i<ns_j; i++){
		if(sizeof(REAL)==sizeof(float)) {
			fprintf(tempfile,"%f %f %f\n%f %f %f\n",
					h_Ra[i],h_Rb[i],h_Rc[i],h_Ba[i],h_Bb[i],h_Bc[i]);
		} else {
			fprintf(tempfile,"%lf %lf %lf\n%lf %lf %lf\n",
					h_Ra[i],h_Rb[i],h_Rc[i],h_Ba[i],h_Bb[i],h_Bc[i]);
		}
	}
	
	fflush(tempfile);
	
	fclose(tempfile);
	
}
//-------------------------------------------------------------------------------------------------------
// Dump_J
//-------------------------------------------------------------------------------------------------------
void Dump_J(void){
	
	FILE* tempfile;
	
	tempfile = Fopen("j_dump_start","w");
	
	if(sizeof(REAL)==sizeof(float)) {
		
		for(int i = 0; i<ns_j; i++) {
			fprintf(tempfile,"%f %f %f\n",h_Jpx[i],h_Jpy[i],h_Jpz[i]);
			fprintf(tempfile,"%f %f %f\n",h_Jmx[i],h_Jmy[i],h_Jmz[i]);
		}
		
	}else{
		
		for(int i = 0; i<ns_j; i++) {
			fprintf(tempfile,"%lf %lf %lf\n",h_Jpx[i],h_Jpy[i],h_Jpz[i]);
			fprintf(tempfile,"%lf %lf %lf\n",h_Jmx[i],h_Jmy[i],h_Jmz[i]);
		}
		
	}
	
	
	fflush(tempfile);
	
	fclose(tempfile);
	
}


//-------------------------------------------------------------------------------------------------------
// Initialize_Spin
//-------------------------------------------------------------------------------------------------------
void Initialize_Spin(REAL* sa_a, REAL* sa_b, REAL* sa_c){
	
	REAL a, b, c;
	REAL invnorma, norma;
	
	for(int i = 0; i<h_nspin; i++){
		
		norma=two;
		
		while (norma>one){
			
			a=(RANDOM01)*two-one;
			b=(RANDOM01)*two-one;
			c=(RANDOM01)*two-one;
			
			norma=a*a+b*b+c*c;
			
		}
		
		invnorma=one/SQRT(norma);
		
		sa_a[i]=a*invnorma;
		sa_b[i]=b*invnorma;
		sa_c[i]=c*invnorma;
		
	}
	
}


//-------------------------------------------------------------------------------------------------------
// Read_Spin
//-------------------------------------------------------------------------------------------------------
void Read_Spin(char *inputfile,
               REAL* sr_a, REAL* sr_b, REAL* sr_c,
               REAL* sb_a, REAL* sb_b, REAL* sb_c) {
	
	REAL a, b, c;
	int j=0, k=0;
	FILE *fpi;
	fpi=Fopen(inputfile,"r");
	fprintf(stderr,"Reading spin from file %s\n",inputfile);
	
	for(int i = 0; i< V(h_side); i++){
		if(sizeof(REAL)==sizeof(float)) {
			fscanf(fpi,"%f %f %f",&a,&b,&c);
		} else {
			fscanf(fpi,"%lf %lf %lf",&a,&b,&c);
		}
		if(i&1) {
			sb_a[j]=a;
			sb_b[j]=b;
			sb_c[j]=c;
			j++;
		} else {
			sr_a[k]=a;
			sr_b[k]=b;
			sr_c[k]=c;
			k++;
		}
	}
	fclose(fpi);
	
}


//-------------------------------------------------------------------------------------------------------
// CopyDataOnDevice
//-------------------------------------------------------------------------------------------------------
void CopyDataOnDevice(void){
	
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Ra, ms_s_mpi ) );
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Rb, ms_s_mpi ) );
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Rc, ms_s_mpi ) );
	
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Ba, ms_s_mpi ) );
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Bb, ms_s_mpi ) );
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Bc, ms_s_mpi ) );
	
	MY_CUDA_CHECK( cudaMemcpy( d_Ra+h_Lsdb2, h_Ra, ms_s, cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Rb+h_Lsdb2, h_Rb, ms_s, cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Rc+h_Lsdb2, h_Rc, ms_s, cudaMemcpyHostToDevice ) );
	
	MY_CUDA_CHECK( cudaMemcpy( d_Ba+h_Lsdb2, h_Ba, ms_s, cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Bb+h_Lsdb2, h_Bb, ms_s, cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Bc+h_Lsdb2, h_Bc, ms_s, cudaMemcpyHostToDevice ) );
	
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jpx, ms_j_mpi ) );
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jpy, ms_j_mpi ) );
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jpz, ms_j_mpi ) );
	
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jmx, ms_j_mpi ) );
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jmy, ms_j_mpi ) );
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jmz, ms_j_mpi ) );
	
	MY_CUDA_CHECK( cudaMemcpy( d_Jpx+h_Lsdb2, h_Jpx, ms_j, cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jpy+h_Lsdb2, h_Jpy, ms_j, cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jpz+h_Lsdb2, h_Jpz, ms_j, cudaMemcpyHostToDevice ) );
	
	MY_CUDA_CHECK( cudaMemcpy( d_Jmx+h_Lsdb2, h_Jmx, ms_j, cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jmy+h_Lsdb2, h_Jmy, ms_j, cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jmz+h_Lsdb2, h_Jmz, ms_j, cudaMemcpyHostToDevice ) );
	
	MY_CUDA_CHECK( cudaMalloc( (void **) &d_Energy, ms_s ) );
	MY_CUDA_CHECK( cudaMemset( d_Energy, 0, ms_s ) );
}


//-------------------------------------------------------------------------------------------------------
// Initialize_J
//-------------------------------------------------------------------------------------------------------
void Initialize_J(void){
	
	int j;
	
	switch (init_type_j){
		
		case 0:
		for (j=0;j<ns_j;j++) {
			h_Jpx[j] = ((RANDOM01)>zerofive) ? -1 : 1;
			h_Jpy[j] = ((RANDOM01)>zerofive) ? -1 : 1;
			h_Jpz[j] = ((RANDOM01)>zerofive) ? -1 : 1;
			h_Jmx[j] = ((RANDOM01)>zerofive) ? -1 : 1;
			h_Jmy[j] = ((RANDOM01)>zerofive) ? -1 : 1;
			h_Jmz[j] = ((RANDOM01)>zerofive) ? -1 : 1;
		}
		
		break;
		case 1:
		for (j=0;j<ns_j;j++) {
			h_Jpx[j] = h_Jpy[j] = h_Jpz[j] =
			h_Jmx[j] = h_Jmy[j] = h_Jmz[j] = one;
		}
		break;
		case -1:
		for (j=0;j<ns_j;j++) {
			h_Jpx[j] = h_Jpy[j] = h_Jpz[j] =
			h_Jmx[j] = h_Jmy[j] = h_Jmz[j] = -one;
		}
		
		break;
		case 2: //Caso Gaussiano
		for (j=0;j<ns_j;j++) {
			h_Jpx[j] = rand_gauss();
			h_Jpy[j] = rand_gauss();
			h_Jpz[j] = rand_gauss();
			h_Jmx[j] = rand_gauss();
			h_Jmy[j] = rand_gauss();
			h_Jmz[j] = rand_gauss();
		}
		break;
		default: {
			writelog(TRUE,APPLICATION_RC,"Invalid initialization for J: %d", init_type_j);
			break;
		}
	}
	
}


//-------------------------------------------------------------------------------------------------------
// rand_gauss
//-------------------------------------------------------------------------------------------------------
REAL rand_gauss(void){
	
	REAL r,cosr;
	REAL a=one-RANDOM01;
	if(a<EPSILON) { a=EPSILON; }
	
	r=SQRT(-2.*LOG(a));
	
	cosr=COS((2*M_PI)*(RANDOM01));
	
	return r*cosr;
	
}


//-------------------------------------------------------------------------------------------------------
// Parameters_Reading
//-------------------------------------------------------------------------------------------------------
void Parameters_Reading(char *inputfile){
	
	dictionary *ini;
	char key[MAXSTRLEN];
	
	if(inputfile==NULL) { writelog(TRUE,APPLICATION_RC,"No inputfile"); }
	
	ini = iniparser_load(inputfile);
	
	if(ini==NULL) { writelog(TRUE,APPLICATION_RC,"Cannot parse file: %s\n", inputfile); }
	
	READINTFI(h_side,"Side");
	READINTFI(init_type_j,"Initialization type for J");
	READINTFI(nsteps,"Number of steps");
	READINTFI(overrelaxiter,"Number of overrelaxation iterations");
	READINTFI(outfreq,"Output frequency");
	READINTFI(dumpfreq,"Dump configuration frequency");
	READREALFI(h_h,"Magnetic field");
	READREALFI(h_Beta,"Beta (1/T)");
	READINTFI(seed,"Random seed");
	
	DEVICEREADINTFI(numBlocks,"Number of blocks");
	DEVICEREADINTFI(numThreadsPerBlock,"Number of threads for block");
	DEVICEREADINTFI(whichgpu,"Which GPU");
	
	DEVICEREADINTFI(RandNumBlocks,"RandGene num of blocks");
	DEVICEREADINTFI(RandNumThreadsPerBlock,"RandGene num of threads for block");
	
}


//-------------------------------------------------------------------------------------------------------
// Calculate_Sizes
//-------------------------------------------------------------------------------------------------------
void Calculate_Sizes(void){
	
	ns_j = h_side*h_side*z_side/2;
	
	ms_s = (h_side*h_side*z_side/2)*sizeof(REAL);
	ms_s_mpi = ((h_side*h_side*z_side/2)+(h_side*h_side))*sizeof(REAL);
	
	ms_j = ns_j*sizeof(REAL);
	ms_j_mpi = ((h_side*h_side*z_side/2)+(h_side*h_side))*sizeof(REAL);
	
}


//-------------------------------------------------------------------------------------------------------
// getNumBlocksAndThreads
//-------------------------------------------------------------------------------------------------------
void getNumBlocksAndThreads(int n, int maxBlocks, int maxThreads, int &blocks, int &threads) {
	
	if (n == 1){
		threads = 1;
		blocks = 1;
	} else{
		threads = (n < maxThreads*2) ? nextPow2(n / 2) : maxThreads;
		blocks = max(1, n / (threads * 2));
	}
	
	blocks = min(maxBlocks, blocks);
}


//-------------------------------------------------------------------------------------------------------
// main
//-------------------------------------------------------------------------------------------------------
int main(int argc, char **argv){
	
	char *inputfile = NULL;
	char *po;
	int istep;
	int nblocks_red, nthreads_red;
	int cudadevice, ngpu;
	double ttimeover=0., ttimehb=0.;
	REAL normene;
	int j;
	char *readspin=NULL;
#if !defined(HOSTREDU)
	unsigned int retCnt;
#endif
	int myid=0, nprocs=1;
	int lpoff;
	struct mpiexch *dataexch;
	int numBlocksBnd;
	
	cudaStream_t *stream;
	int nstreams=2;
	
	TIMER_DEF;
	
	for(int i = 1; i < argc; i++) {
		po = argv[i];
		if (*po++ == '-') {
			switch (*po++) {
				case '?':
				SKIPBLANK
				Usage(argv[0]);
				exit(OK);
				break;
				case 'v':
				SKIPBLANK
				verbose=TRUE;
				break;
				case 'i':
				SKIPBLANK
				inputfile=Strdup(po);
				break;
				case 'o':
				SKIPBLANK
				onlyOverrelaxation=true;
				break;
				case 'r':
				SKIPBLANK
				readspin=Strdup(po);
				break;
				
				case 'b':
				SKIPBLANK
				onlyHeatBath=true;
				break;
				default:
				Usage(argv[0]);
				exit(OK);
				break;
			}
		}
	}
	if(inputfile==NULL) {
		Usage(argv[0]);
		exit(OK);
	}
	
	StartMpi(&myid, &nprocs,&argc,&argv);
	mpiid=myid;
	
#if defined(DEBUG)
	cudaPrintfInit();
#endif

	cudaGetDeviceCount(&ngpu);
	
	if(ngpu<1) {
		writelog(TRUE,APPLICATION_RC,"No gpu available! %d\n",ngpu);
	}
	
	Parameters_Reading(inputfile);
	
	seed=(seed+myid*nprocs)%MAXSEED;
	
	
	z_side=h_side/nprocs;
	
	numBlocksBnd=h_Lsdb2/numThreadsPerBlock+(h_Lsdb2%numThreadsPerBlock?1:0);
	
	if((z_side*nprocs)<h_side) {
		writelog(TRUE,APPLICATION_RC,"Invalid number of processors (%d) or system size (%d). Number of processors must be an exact divider of the system size", nprocs, h_side);
	}
	
	if(whichgpu>(ngpu-1)) {
        	writelog(TRUE,APPLICATION_RC, 
		"invalid number of gpu: %d must be > 0 and <= %d\n",ngpu);
	}
	
	
	if(nprocs>1) {
		if(myid%2) { /* Odd task go on GPU 1 */
			whichgpu=1;
		} else {
			whichgpu=0;
		}
	}

	if(whichgpu>=0) {
        	printf("There are %d GPUs!, forcing usage of n. <%d>\n",ngpu,whichgpu);
        	MY_CUDA_CHECK(cudaSetDevice(whichgpu));
	}

	dataexch=SetUpMpi(myid,nprocs,h_side);
	
	if (h_side<=0 || !pow_of_two(h_side)) {
		writelog(TRUE,APPLICATION_RC,"Invalid side: must be > 0 and power of two\n");
	}
	
	Calculate_Sizes();
	Init();
	
	lpoff=(z_side-1)*h_Lsdb2;
	
	if(readspin==NULL) {
		
		Initialize_Spin(h_Ra,h_Rb,h_Rc);
		Initialize_Spin(h_Ba,h_Bb,h_Bc);
		Dump_Spin(-1);
		
	} else {
		Read_Spin(readspin,h_Ra,h_Rb,h_Rc, h_Ba,h_Bb,h_Bc);
	}
	
	Initialize_J();
	
	
	// Dump_J();
	
	initRand(V(h_side));
	
	CopyConstOnDevice();
	
	CopyDataOnDevice();
	
	getNumBlocksAndThreads(h_nspin, numBlocks, numThreadsPerBlock,nblocks_red, nthreads_red);
	
	normene=one/(three*V(h_side));
	
	h_totene=zeroEne;
	
	if(myid==0 && nprocs==1) {
		ComputeHEnergy(h_Ra, h_Rb, h_Rc, h_Ba, h_Bb, h_Bc,h_Jpx, h_Jpy, h_Jpz, h_Jmx, h_Jmy, h_Jmz, &h_totene,TRUE);
		h_totene=SumEnergy(h_totene, myid);
		printf("Initial Energy per spin on host: %f\n", normene*h_totene);
	}
	
	memcpy(dataexch->downsend,h_Ra,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->downsend+h_Lsdb2,h_Rb,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->downsend+(2*h_Lsdb2),h_Rc,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->downsend+(3*h_Lsdb2),h_Ba,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->downsend+(4*h_Lsdb2),h_Bb,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->downsend+(5*h_Lsdb2),h_Bc,h_Lsdb2*sizeof(REAL));
	
	memcpy(dataexch->upsend,h_Ra+lpoff,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->upsend+h_Lsdb2,h_Rb+lpoff,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->upsend+(2*h_Lsdb2),h_Rc+lpoff,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->upsend+(3*h_Lsdb2),h_Ba+lpoff,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->upsend+(4*h_Lsdb2),h_Bb+lpoff,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->upsend+(5*h_Lsdb2),h_Bc+lpoff,h_Lsdb2*sizeof(REAL));
	
	
	dataexch->nbyte=2*h_Lsdb2*sizeof(REAL)*NCOMP;
	dataexch->post=FALSE;
	
	ExchMpi(dataexch);
	
	MY_CUDA_CHECK( cudaMemcpy( d_Ra, dataexch->downrecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Rb, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Rc, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Ba, dataexch->downrecv+(3*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Bb, dataexch->downrecv+(4*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Bc, dataexch->downrecv+(5*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	
	MY_CUDA_CHECK( cudaMemcpy( d_Ra+h_nspin+h_Lsdb2, dataexch->uprecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Rb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Rc+h_nspin+h_Lsdb2, dataexch->uprecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Ba+h_nspin+h_Lsdb2, dataexch->uprecv+(3*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Bb+h_nspin+h_Lsdb2, dataexch->uprecv+(4*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Bc+h_nspin+h_Lsdb2, dataexch->uprecv+(5*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	
	memcpy(dataexch->downsend,h_Jpx,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->downsend+h_Lsdb2,h_Jpy,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->downsend+(2*h_Lsdb2),h_Jpz,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->downsend+(3*h_Lsdb2),h_Jmx,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->downsend+(4*h_Lsdb2),h_Jmy,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->downsend+(5*h_Lsdb2),h_Jmz,h_Lsdb2*sizeof(REAL));
	
	memcpy(dataexch->upsend,h_Jpx+lpoff,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->upsend+h_Lsdb2,h_Jpy+lpoff,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->upsend+(2*h_Lsdb2),h_Jpz+lpoff,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->upsend+(3*h_Lsdb2),h_Jmx+lpoff,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->upsend+(4*h_Lsdb2),h_Jmy+lpoff,h_Lsdb2*sizeof(REAL));
	memcpy(dataexch->upsend+(5*h_Lsdb2),h_Jmz+lpoff,h_Lsdb2*sizeof(REAL));
	
	dataexch->post=FALSE;
	
	ExchMpi(dataexch);
	
	MY_CUDA_CHECK( cudaMemcpy( d_Jpx, dataexch->downrecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jpy, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jpz, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jmx, dataexch->downrecv+(3*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jmy, dataexch->downrecv+(4*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jmz, dataexch->downrecv+(5*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	
	MY_CUDA_CHECK( cudaMemcpy( d_Jpx+h_nspin+h_Lsdb2, dataexch->uprecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jpy+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jpz+h_nspin+h_Lsdb2, dataexch->uprecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jmx+h_nspin+h_Lsdb2, dataexch->uprecv+(3*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jmy+h_nspin+h_Lsdb2, dataexch->uprecv+(4*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	MY_CUDA_CHECK( cudaMemcpy( d_Jmz+h_nspin+h_Lsdb2, dataexch->uprecv+(5*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	
	dataexch->post=FALSE;
	
	dataexch->nbyte=h_Lsdb2*sizeof(REAL)*NCOMP;
	
	MY_CUDA_CHECK( cudaMalloc((void **)&r_odata,numBlocks*sizeof(REAL)) );
	
	dim3 dimGrid(numBlocks,1,1);
	dim3 dimGridBnd(numBlocksBnd,1,1);
	dim3 dimBlock(numThreadsPerBlock,1,1);
	
	ComputeEnergy<<< dimGrid,dimBlock >>>(d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
										  d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
										  d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
										  d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
										  d_Energy, TRUE);
										  
#if defined(HOSTREDU)
	MY_CUDA_CHECK( cudaMemcpy(h_Energy, d_Energy, h_nspin*sizeof(REAL),cudaMemcpyDeviceToHost) );
	
	h_totene=zeroEne;
	
	for(int i=0; i<h_nspin; i++) { h_totene+=h_Energy[i]; }
	
#else

	MY_CUDA_CHECK( cudaThreadSynchronize() );
	retCnt=0;
	cudaMemcpyToSymbol("retirementCount", &retCnt, sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
	reduceSinglePass(h_nspin, nthreads_red, nblocks_red, d_Energy, r_odata);
	MY_CUDA_CHECK( cudaThreadSynchronize() );
	MY_CUDA_CHECK( cudaMemcpy(&h_totene, r_odata, sizeof(REAL),cudaMemcpyDeviceToHost) );
	
#endif

	h_totene=SumEnergy(h_totene, myid);
	if(myid==0) {
		printf("Initial Energy per spin: %f\n", normene*h_totene);
	}
	
#if defined(FERMI)
	cudaFuncSetCacheConfig(OverrelaxRedBnd, cudaFuncCachePreferL1 );
	cudaFuncSetCacheConfig(OverrelaxRedBulk, cudaFuncCachePreferL1 );
	cudaFuncSetCacheConfig(OverrelaxBlueBnd, cudaFuncCachePreferL1 );
	cudaFuncSetCacheConfig(OverrelaxBlueBulk, cudaFuncCachePreferL1 );
#endif
	// allocate and initialize an array of stream handles
	stream = (cudaStream_t*) Malloc(nstreams * sizeof(cudaStream_t));
	
	for(int i = 0; i < nstreams; i++) {
		MY_CUDA_CHECK( cudaStreamCreate(&(stream[i])) );
	}
		
	for(istep=0; istep<nsteps; istep++) {
		
		if(!onlyHeatBath){
			
			TIMER_START;
			
			for(j=0; j<overrelaxiter; j++) {
				
#if defined(SINGLEOVERKERNEL)

				Overrelaxation<<< dimGrid,dimBlock>>>(d_Ra, d_Rb, d_Rc,
													  d_Ba, d_Bb, d_Bc,
													  d_Jpx, d_Jpy, d_Jpz,
													  d_Jmx, d_Jmy, d_Jmz, TRUE);
				
				MY_CUDA_CHECK( cudaThreadSynchronize() );
				
				Overrelaxation<<< dimGrid,dimBlock>>>(d_Ba, d_Bb, d_Bc,
													  d_Ra, d_Rb, d_Rc,
													  d_Jmx, d_Jmy, d_Jmz,
													  d_Jpx, d_Jpy, d_Jpz, FALSE);
				
				MY_CUDA_CHECK( cudaThreadSynchronize() );
				
#else

				OverrelaxRedBnd<<< dimGridBnd,dimBlock, 0, stream[1] >>>(d_Ra+h_Lsdb2,
																		 d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
																		 d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
																		 d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
																		 d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2);
				OverrelaxRedBulk<<< dimGrid,dimBlock, 0, stream[0] >>>(d_Ra+h_Lsdb2,
																	   d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
																	   d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
																	   d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
																	   d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2);
																	   
				MY_CUDA_CHECK(cudaStreamSynchronize(stream[0]));
				
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend, d_Ra+h_Lsdb2,h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+h_Lsdb2, d_Rb+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+(2*h_Lsdb2), d_Rc+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend, d_Ra+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+h_Lsdb2, d_Rb+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+(2*h_Lsdb2), d_Rc+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
				MY_CUDA_CHECK(cudaStreamSynchronize(stream[1]));
				
				ExchMpi(dataexch);
				
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Ra, dataexch->downrecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Rb, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Rc, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Ra+h_nspin+h_Lsdb2, dataexch->uprecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Rb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Rc+h_nspin+h_Lsdb2,dataexch->uprecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				
				MY_CUDA_CHECK( cudaThreadSynchronize() );
				
				OverrelaxBlueBnd<<< dimGridBnd,dimBlock, 0, stream[1]>>>(d_Ba+h_Lsdb2,
																		 d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
																		 d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
																		 d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
																		 d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2);
				OverrelaxBlueBulk<<< dimGrid,dimBlock, 0, stream[0]>>>(d_Ba+h_Lsdb2,
																	   d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
																	   d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
																	   d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
																	   d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2);
																	   
				MY_CUDA_CHECK(cudaStreamSynchronize(stream[0]));
				
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend, d_Ba+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+h_Lsdb2, d_Bb+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+(2*h_Lsdb2), d_Bc+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend, d_Ba+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+h_Lsdb2, d_Bb+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+(2*h_Lsdb2), d_Bc+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
				MY_CUDA_CHECK(cudaStreamSynchronize(stream[1]));
				
				ExchMpi(dataexch);
				
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Ba, dataexch->downrecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Bb, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Bc, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Ba+h_nspin+h_Lsdb2, dataexch->uprecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Bb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Bc+h_nspin+h_Lsdb2,dataexch->uprecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaThreadSynchronize() );
				
#endif
			}
			
			TIMER_STOP;
			ttimeover+=TIMER_ELAPSED;
		}
		
		if(!onlyOverrelaxation) {
			
			TIMER_START;
			
			/* Here it goes the heat-bath */
			/* first of all, random number generation */
			/* then the heat-bath itself */
			ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ONEONE,d_cRand);
			ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ZEROFIVE,d_aRand);
			ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ZEROFIVE,d_bRand);
			ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ZEROONE,d_Rand);
			      
			if(istep==0&&FALSE){
				printRand(randNum,d_aRand,"daRand3");
				printRand(randNum,d_bRand,"dbRand3");
				printRand(randNum,d_cRand,"dcRand3");
				printRand(randNum,d_Rand,"dRand3");
				exit(1);
			}
			
			
#if defined(SINGLEHBKERNEL)

			HeatBath<<< dimGrid,dimBlock >>>(d_Ra, d_Rb, d_Rc,
											 d_Ba, d_Bb, d_Bc,
											 d_Jpx, d_Jpy, d_Jpz,
											 d_Jmx, d_Jmy, d_Jmz,
											 d_aRand,d_bRand,d_cRand,
											 d_Rand,TRUE);
			
			MY_CUDA_CHECK( cudaThreadSynchronize() );
			
			HeatBath<<< dimGrid,dimBlock >>>(d_Ba, d_Bb, d_Bc,
											 d_Ra, d_Rb, d_Rc,
											 d_Jmx, d_Jmy, d_Jmz,
											 d_Jpx, d_Jpy, d_Jpz,
											 d_aRand+h_nspin,d_bRand+h_nspin,
											 d_cRand+h_nspin,
											 d_Rand+h_nspin,FALSE);
			
			MY_CUDA_CHECK( cudaThreadSynchronize() );
      
#else


				HeatBathRedBnd<<< dimGridBnd,dimBlock, 0, stream[1] >>>(d_Ra+h_Lsdb2,
																		 d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
																		 d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
																		 d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
																		 d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
																		d_aRand,d_bRand,d_cRand,d_Rand);
							
				HeatBathRedBulk<<< dimGrid,dimBlock, 0, stream[0] >>>(d_Ra+h_Lsdb2,
																	   d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
																	   d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
																	   d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
																	   d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
																		d_aRand,d_bRand,d_cRand,d_Rand);
																	   
				MY_CUDA_CHECK(cudaStreamSynchronize(stream[0]));

				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend, d_Ra+h_Lsdb2,h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+h_Lsdb2, d_Rb+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+(2*h_Lsdb2), d_Rc+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend, d_Ra+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+h_Lsdb2, d_Rb+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+(2*h_Lsdb2), d_Rc+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
				MY_CUDA_CHECK(cudaStreamSynchronize(stream[1]));
				
				ExchMpi(dataexch);
				
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Ra, dataexch->downrecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Rb, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Rc, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Ra+h_nspin+h_Lsdb2, dataexch->uprecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Rb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Rc+h_nspin+h_Lsdb2,dataexch->uprecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				
				MY_CUDA_CHECK( cudaThreadSynchronize() );
				
				HeatBathBlueBnd<<< dimGridBnd,dimBlock, 0, stream[1]>>>(d_Ba+h_Lsdb2,
																		 d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
																		 d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
																		 d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
																		 d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
																		d_aRand,d_bRand,d_cRand,d_Rand);

				HeatBathBlueBulk<<< dimGrid,dimBlock, 0, stream[0]>>>(d_Ba+h_Lsdb2,
																	   d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
																	   d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
																	   d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
																	   d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
																		d_aRand,d_bRand,d_cRand,d_Rand);
																	   
				MY_CUDA_CHECK(cudaStreamSynchronize(stream[0]));

				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend, d_Ba+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+h_Lsdb2, d_Bb+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+(2*h_Lsdb2), d_Bc+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend, d_Ba+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+h_Lsdb2, d_Bb+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+(2*h_Lsdb2), d_Bc+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
				MY_CUDA_CHECK(cudaStreamSynchronize(stream[1]));
				
				ExchMpi(dataexch);
				
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Ba, dataexch->downrecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Bb, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Bc, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Ba+h_nspin+h_Lsdb2, dataexch->uprecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Bb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				MY_CUDA_CHECK( cudaMemcpyAsync( d_Bc+h_nspin+h_Lsdb2,dataexch->uprecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );

				MY_CUDA_CHECK( cudaThreadSynchronize() );




	/*		HeatBathRed<<< dimGrid,dimBlock >>>(d_Ra, d_Rb, d_Rc,
												d_Ba, d_Bb, d_Bc,
												d_Jpx, d_Jpy, d_Jpz,
												d_Jmx, d_Jmy, d_Jmz,
												d_aRand,d_bRand,d_cRand,
												d_Rand);
			
			MY_CUDA_CHECK( cudaThreadSynchronize() );
			
			HeatBathBlue<<< dimGrid,dimBlock >>>(d_Ba, d_Bb, d_Bc,
												 d_Ra, d_Rb, d_Rc,
												 d_Jmx, d_Jmy, d_Jmz,
												 d_Jpx, d_Jpy, d_Jpz,
												 d_aRand+h_nspin,d_bRand+h_nspin,
												 d_cRand+h_nspin,
												 d_Rand+h_nspin);
			
			MY_CUDA_CHECK( cudaThreadSynchronize() );*/
			
#endif
			
			TIMER_STOP;
			ttimehb+=TIMER_ELAPSED;
			
		}
		
		if(istep && (istep%outfreq)==0) {
			ComputeEnergy<<< dimGrid,dimBlock >>>(d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
												  d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
												  d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
												  d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
												  d_Energy, TRUE);
			
			MY_CUDA_CHECK( cudaThreadSynchronize() );
			
			
#if defined(HOSTREDU)
			
			MY_CUDA_CHECK( cudaMemcpy(h_Energy, d_Energy,
									  h_nspin*sizeof(REAL),
									  cudaMemcpyDeviceToHost) );
			h_totene=zeroEne;
			
			for(int i=0; i<h_nspin; i++) h_totene+=h_Energy[i];
			
#else
			retCnt=0;
			cudaMemcpyToSymbol("retirementCount", &retCnt, sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
			reduceSinglePass(h_nspin, nthreads_red, nblocks_red, d_Energy, r_odata);
			MY_CUDA_CHECK( cudaMemcpy(&h_totene, r_odata, sizeof(REAL),cudaMemcpyDeviceToHost) );
			
#endif
			h_totene=SumEnergy(h_totene, myid);
			if(myid==0) {
				printf("Energy per spin at step %d: %f\n",istep,normene*h_totene);
			}
			//place holder for energy evaluation
		}
		
		if(istep && (istep%dumpfreq)==0) {
			//Dump_Spin(istep);
		}
		
	}
	
	if(myid==0) {
		if(!onlyHeatBath){
			printf("Total time for %d iterations of overrelaxation on a %d^3 cube: %f seconds\n", nsteps*overrelaxiter, h_side, ttimeover*0.000001);
			printf("Time for spin: %f microseconds\n",ttimeover/(overrelaxiter*(double)nsteps*V(h_side)));
		}
		
		if(!onlyOverrelaxation){
			printf("Total time for %d iterations of heat bath on a %d^3 cube: %f seconds\n", nsteps, h_side, ttimehb*0.000001);
			printf("Time for spin: %f microseconds\n",ttimehb/((double)nsteps*V(h_side)));
		}
	}
	
	ComputeEnergy<<< dimGrid,dimBlock >>>(d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
										  d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
										  d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
										  d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
										  d_Energy, TRUE);
	
#if defined(HOSTREDU)
	
	MY_CUDA_CHECK( cudaMemcpy(h_Energy, d_Energy, h_nspin*sizeof(REAL),cudaMemcpyDeviceToHost) );
	
	h_totene=zero;
	
	for(int i=0; i<h_nspin; i++) h_totene+=h_Energy[i];
	
#else
	
	MY_CUDA_CHECK( cudaThreadSynchronize() );
	reduceSinglePass(ns_j, nthreads_red, nblocks_red,d_Energy, r_odata);
	MY_CUDA_CHECK( cudaMemcpy(&h_totene, r_odata, sizeof(REAL),cudaMemcpyDeviceToHost) );
	
#endif

	h_totene=SumEnergy(h_totene, myid);
	
	if(myid==0) {
		printf("Final Energy per spin: %f\n",normene*h_totene);
	}
	
	Dump_Spin(nsteps);
	
	FreeDevice();
	
	FreeHost();
	
	FreeRand();
	
#if defined(DEBUG)
	cudaPrintfDisplay();
	cudaPrintfEnd();
#endif

	StopMpi();
	
	return 0;
	
}


//-------------------------------------------------------------------------------------------------------
// Overrelaxation
//-------------------------------------------------------------------------------------------------------
#if defined(SINGLEOVERKERNEL)
__global__ void Overrelaxation(REAL *a, REAL *b, REAL *c,
                               REAL *an, REAL *bn, REAL *cn,
                               REAL *Jpx, REAL *Jpy, REAL *Jpz,
                               REAL *Jmx, REAL *Jmy, REAL *Jmz,
                               int even) {
	
	unsigned int iy, iz;
	int diff;
	REAL as, bs, cs;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart;
	
	istart=blockDim.x*blockIdx.x+d_Lsdb2; /*start on second plane (by using d_Lsdb2)
										   since first plane is for the boundary */
	
	
	for(i=istart+tid; i<d_nspin; i+=num_threads) {
		
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		//do the same thing if both $iz and $iy are even or odd
		diff=(((iz&1)^(iy&1))^even);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2) & d_maskL3;
		smz=(i-d_Lsdb2) & d_maskL3;
		
		as=an[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
		an[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
		an[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
		an[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
		an[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
		an[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];
		
		bs=bn[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
		bn[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
		bn[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
		bn[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
		bn[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
		bn[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];
		
		cs=cn[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
		cn[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
		cn[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
		cn[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
		cn[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
		cn[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];
		
		cs+=d_h;
		
		//Overrelaxation
		
		factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/(as*as+bs*bs+cs*cs);
		
		a[i]=as*factor-a[i];
		b[i]=bs*factor-b[i];
		c[i]=cs*factor-c[i];
		
	}
	
}

#else

//-------------------------------------------------------------------------------------------------------
// OverrelaxRedBnd
//-------------------------------------------------------------------------------------------------------
__global__ void OverrelaxRedBnd(REAL *a, REAL *b, REAL *c,
								REAL *an, REAL *bn, REAL *cn,
								REAL *Jpx, REAL *Jpy, REAL *Jpz,
								REAL *Jmx, REAL *Jmy, REAL *Jmz) {
	
	unsigned int iy, iz;
	int diff;
	REAL as, bs, cs;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart;
	
	istart=blockDim.x*blockIdx.x; /* does first plane (up to d_Lsdb2) */
	
	for(i=istart+tid; i<d_Lsdb2; i+=num_threads) {
		
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		//do the same thing if both $iz and $iy are even or odd
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2) /* & d_maskL3 */;
		smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
		as=an[spx]*Jpx[i]+
		an[spy]*Jpy[i]+
		an[spz]*Jpz[i]+
		an[smx]*Jmx[i]+
		an[smy]*Jmy[i]+
		an[smz]*Jmz[i];
		
		bs=bn[spx]*Jpx[i]+
		bn[spy]*Jpy[i]+
		bn[spz]*Jpz[i]+
		bn[smx]*Jmx[i]+
		bn[smy]*Jmy[i]+
		bn[smz]*Jmz[i];
		
		cs=cn[spx]*Jpx[i]+
		cn[spy]*Jpy[i]+
		cn[spz]*Jpz[i]+
		cn[smx]*Jmx[i]+
		cn[smy]*Jmy[i]+
		cn[smz]*Jmz[i];
		
		cs+=d_h;
		
		//Overrelaxation
		
		factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/(as*as+bs*bs+cs*cs);
		
#if defined(DEBUG)
		printf("Reds [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif
		
		a[i]=as*factor-a[i];
		b[i]=bs*factor-b[i];
		c[i]=cs*factor-c[i];
		
#if defined(DEBUG)
		printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif
		
		
	}
	istart=(d_nspin-d_Lsdb2)+blockDim.x*blockIdx.x;
	
	for(i=istart+tid; i<d_nspin; i+=num_threads) {
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		//do the same thing if both $iz and $iy are even or odd
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2) /* & d_maskL3 */;
		smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
		as=an[spx]*Jpx[i]+
		an[spy]*Jpy[i]+
		an[spz]*Jpz[i]+
		an[smx]*Jmx[i]+
		an[smy]*Jmy[i]+
		an[smz]*Jmz[i];
		
		bs=bn[spx]*Jpx[i]+
		bn[spy]*Jpy[i]+
		bn[spz]*Jpz[i]+
		bn[smx]*Jmx[i]+
		bn[smy]*Jmy[i]+
		bn[smz]*Jmz[i];
		
		cs=cn[spx]*Jpx[i]+
		cn[spy]*Jpy[i]+
		cn[spz]*Jpz[i]+
		cn[smx]*Jmx[i]+
		cn[smy]*Jmy[i]+
		cn[smz]*Jmz[i];
		
		cs+=d_h;
		
		//Overrelaxation
		
		factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/(as*as+bs*bs+cs*cs);
		
#if defined(DEBUG)
		printf("Reds [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif
		
		a[i]=as*factor-a[i];
		b[i]=bs*factor-b[i];
		c[i]=cs*factor-c[i];
		
#if defined(DEBUG)
		printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif
		
	}
}


//-------------------------------------------------------------------------------------------------------
// OverrelaxRedBulk
//-------------------------------------------------------------------------------------------------------
__global__ void OverrelaxRedBulk(REAL *a, REAL *b, REAL *c,
								 REAL *an, REAL *bn, REAL *cn,
								 REAL *Jpx, REAL *Jpy, REAL *Jpz,
								 REAL *Jmx, REAL *Jmy, REAL *Jmz) {
	
	unsigned int iy, iz;
	int diff;
	REAL as, bs, cs;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart, iend;
	
	istart=blockDim.x*blockIdx.x+d_Lsdb2; /*start on second plane (by using d_Lsdb2)
										   since first plane is for the boundary */
	iend=d_nspin-d_Lsdb2;
	
	for(i=istart+tid; i<iend; i+=num_threads) {
		
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		//do the same thing if both $iz and $iy are even or odd
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2) /* & d_maskL3 */;
		smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
		as=an[spx]*Jpx[i]+
		an[spy]*Jpy[i]+
		an[spz]*Jpz[i]+
		an[smx]*Jmx[i]+
		an[smy]*Jmy[i]+
		an[smz]*Jmz[i];
		
		bs=bn[spx]*Jpx[i]+
		bn[spy]*Jpy[i]+
		bn[spz]*Jpz[i]+
		bn[smx]*Jmx[i]+
		bn[smy]*Jmy[i]+
		bn[smz]*Jmz[i];
		
		cs=cn[spx]*Jpx[i]+
		cn[spy]*Jpy[i]+
		cn[spz]*Jpz[i]+
		cn[smx]*Jmx[i]+
		cn[smy]*Jmy[i]+
		cn[smz]*Jmz[i];
		
		cs+=d_h;
		
		//Overrelaxation
		
		factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/(as*as+bs*bs+cs*cs);
		
#if defined(DEBUG)
		printf("Reds [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif
		
		a[i]=as*factor-a[i];
		b[i]=bs*factor-b[i];
		c[i]=cs*factor-c[i];
		
#if defined(DEBUG)
		printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif
		
		
	}
	
}


//-------------------------------------------------------------------------------------------------------
// OverrelaxBlueBnd
//-------------------------------------------------------------------------------------------------------
__global__ void OverrelaxBlueBnd(REAL *a, REAL *b, REAL *c,
								 REAL *an, REAL *bn, REAL *cn,
								 REAL *Jpx, REAL *Jpy, REAL *Jpz,
								 REAL *Jmx, REAL *Jmy, REAL *Jmz) {
	
	unsigned int iy, iz;
	int diff;
	REAL as, bs, cs;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart;
	
	istart=blockDim.x*blockIdx.x; /* does first plane (by using d_Lsdb2) */
	
	for(i=istart+tid; i<d_Lsdb2; i+=num_threads) {
		
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		//do the same thing if both $iz and $iy are even or odd
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2) /* & d_maskL3 */;
		smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
		as=an[spx]*Jpx[spx]+
		an[spy]*Jpy[spy]+
		an[spz]*Jpz[spz]+
		an[smx]*Jmx[smx]+
		an[smy]*Jmy[smy]+
		an[smz]*Jmz[smz];
		
		bs=bn[spx]*Jpx[spx]+
		bn[spy]*Jpy[spy]+
		bn[spz]*Jpz[spz]+
		bn[smx]*Jmx[smx]+
		bn[smy]*Jmy[smy]+
		bn[smz]*Jmz[smz];
		
		cs=cn[spx]*Jpx[spx]+
		cn[spy]*Jpy[spy]+
		cn[spz]*Jpz[spz]+
		cn[smx]*Jmx[smx]+
		cn[smy]*Jmy[smy]+
		cn[smz]*Jmz[smz];
		
		cs+=d_h;
		
		//Overrelaxation
		
		factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/
		(as*as+bs*bs+cs*cs);
#if defined(DEBUG)
		printf("Blues [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif
		a[i]=as*factor-a[i];
		b[i]=bs*factor-b[i];
		c[i]=cs*factor-c[i];
#if defined(DEBUG)
		printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif
		
	}
	
	istart=(d_nspin-d_Lsdb2)+blockDim.x*blockIdx.x;
	
	for(i=istart+tid; i<d_nspin; i+=num_threads) {
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		//do the same thing if both $iz and $iy are even or odd
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2) /* & d_maskL3 */;
		smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
		as=an[spx]*Jpx[spx]+
		an[spy]*Jpy[spy]+
		an[spz]*Jpz[spz]+
		an[smx]*Jmx[smx]+
		an[smy]*Jmy[smy]+
		an[smz]*Jmz[smz];
		
		bs=bn[spx]*Jpx[spx]+
		bn[spy]*Jpy[spy]+
		bn[spz]*Jpz[spz]+
		bn[smx]*Jmx[smx]+
		bn[smy]*Jmy[smy]+
		bn[smz]*Jmz[smz];
		
		cs=cn[spx]*Jpx[spx]+
		cn[spy]*Jpy[spy]+
		cn[spz]*Jpz[spz]+
		cn[smx]*Jmx[smx]+
		cn[smy]*Jmy[smy]+
		cn[smz]*Jmz[smz];
		
		cs+=d_h;
		
		//Overrelaxation
		
		factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/
		(as*as+bs*bs+cs*cs);
#if defined(DEBUG)
		printf("Blues [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif
		a[i]=as*factor-a[i];
		b[i]=bs*factor-b[i];
		c[i]=cs*factor-c[i];
#if defined(DEBUG)
		printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif
	}
	
}


//-------------------------------------------------------------------------------------------------------
// OverrelaxBlueBulk
//-------------------------------------------------------------------------------------------------------
__global__ void OverrelaxBlueBulk(REAL *a, REAL *b, REAL *c,
								  REAL *an, REAL *bn, REAL *cn,
								  REAL *Jpx, REAL *Jpy, REAL *Jpz,
								  REAL *Jmx, REAL *Jmy, REAL *Jmz) {
	
	unsigned int iy, iz;
	int diff;
	REAL as, bs, cs;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart, iend;
	
	istart=blockDim.x*blockIdx.x; /*start on second plane (by using d_Lsdb2)
								   since first plane is for the boundary */
	iend=d_nspin-d_Lsdb2;
	
	for(i=istart+tid; i<iend; i+=num_threads) {
		
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		//do the same thing if both $iz and $iy are even or odd
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2) /* & d_maskL3 */;
		smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
		as=an[spx]*Jpx[spx]+
		an[spy]*Jpy[spy]+
		an[spz]*Jpz[spz]+
		an[smx]*Jmx[smx]+
		an[smy]*Jmy[smy]+
		an[smz]*Jmz[smz];
		
		bs=bn[spx]*Jpx[spx]+
		bn[spy]*Jpy[spy]+
		bn[spz]*Jpz[spz]+
		bn[smx]*Jmx[smx]+
		bn[smy]*Jmy[smy]+
		bn[smz]*Jmz[smz];
		
		cs=cn[spx]*Jpx[spx]+
		cn[spy]*Jpy[spy]+
		cn[spz]*Jpz[spz]+
		cn[smx]*Jmx[smx]+
		cn[smy]*Jmy[smy]+
		cn[smz]*Jmz[smz];
		
		cs+=d_h;
		
		//Overrelaxation
		
		factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/
		(as*as+bs*bs+cs*cs);
#if defined(DEBUG)
		printf("Blues [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif
		a[i]=as*factor-a[i];
		b[i]=bs*factor-b[i];
		c[i]=cs*factor-c[i];
#if defined(DEBUG)
		printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif
		
	}
	
}


#endif


//-------------------------------------------------------------------------------------------------------
// ComputeEnergy
//-------------------------------------------------------------------------------------------------------
__global__ void ComputeEnergy(REAL *a, REAL *b, REAL *c,
                              REAL *an, REAL *bn, REAL *cn,
                              REAL *Jpx, REAL *Jpy, REAL *Jpz,
                              REAL *Jmx, REAL *Jmy, REAL *Jmz,
                              REAL *energy,
                              int even) {
	
	unsigned int iy, iz;
	int diff;
	int spx, smx, spy, smy, spz, smz;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart;
	
	istart=blockDim.x*blockIdx.x; /*start on second plane (by using d_Lsdb2)
								   since first plane is for the boundary */
	
	
	for(i=istart+tid; i<d_nspin; i+=num_threads) {
		
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		//do the same thing if both $iz and $iy are even or odd
		diff=(((iz&1)^(iy&1))^even);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2) /* & d_maskL3 */;
		smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
		
		energy[i]=a[i]*(an[spx]*Jpx[i]+
						an[spy]*Jpy[i]+
						an[spz]*Jpz[i]+
						an[smx]*Jmx[i]+
						an[smy]*Jmy[i]+
						an[smz]*Jmz[i])+
		b[i]*(bn[spx]*Jpx[i]+
			  bn[spy]*Jpy[i]+
			  bn[spz]*Jpz[i]+
			  bn[smx]*Jmx[i]+
			  bn[smy]*Jmy[i]+
			  bn[smz]*Jmz[i])+
		c[i]*(cn[spx]*Jpx[i]+
			  cn[spy]*Jpy[i]+
			  cn[spz]*Jpz[i]+
			  cn[smx]*Jmx[i]+
			  cn[smy]*Jmy[i]+
			  cn[smz]*Jmz[i]);
#if 0
		cuPrintf("i=%d, spx=%d, spy=%d, spz=%d, smx=%d, smy=%d, smz=%d, ", i, spx, spy, spz, smx, smy, smz);
		cuPrintf("a=%f, b=%f, c=%f, anspx=%f, anspy=%f, anspz=%f, ansmx=%f, ansmy=%f, ansmz=%f, ", a[i], b[i], c[i], an[spx],an[spy],an[spz], an[smx],an[smy],an[smz]);
		
		cuPrintf("bnspx=%f, bnspy=%f, bnspz=%f, bnsmx=%f, bnsmy=%f, bnsmz=%f, ", bn[spx],bn[spy],bn[spz], bn[smx],bn[smy],bn[smz]);
		cuPrintf("cnspx=%f, cnspy=%f, cnspz=%f, cnsmx=%f, cnsmy=%f, cnsmz=%f, ", cn[spx],cn[spy],cn[spz], cn[smx],cn[smy],cn[smz]);
		cuPrintf("Jpx=%f, Jpy=%f, Jpz=%f, Jmx=%f, Jmy=%f, Jmz=%f\n" , Jpx[i],Jpy[i],Jpz[i], Jmx[i],Jmy[smy],Jmz[i]);
#endif
	}
	
}


//-------------------------------------------------------------------------------------------------------
// ComputeHEnergy
//-------------------------------------------------------------------------------------------------------
void ComputeHEnergy(REAL *a, REAL *b, REAL *c,
                    REAL *an, REAL *bn, REAL *cn,
                    REAL *Jpx, REAL *Jpy, REAL *Jpz,
                    REAL *Jmx, REAL *Jmy, REAL *Jmz,
                    REALENE *energy,
                    int even) {
	
	unsigned int iy, iz;
	int diff;
	int spx, smx, spy, smy, spz, smz;
	REAL e;
	
	int i;
	
	for(i=0; i<h_nspin; i++) {
		iz=i>>h_sdbLsm1;
		iy=(i - (iz<<h_sdbLsm1))>>h_sdbLm1;
		
		//do the same thing if both $iz and $iy are even or odd
		diff=(((iz&1)^(iy&1))^even);
		
		spx=(i >> (diff-1)) + (((i-(i&h_maskL)) + ((i+1)&h_maskL)) >> (-diff));
		smx=(((i-(i&h_maskL)) + ((i-1)&h_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&h_maskL2))+((i+h_Ldb2)&h_maskL2);
		smy=(i-(i&h_maskL2))+((i-h_Ldb2)&h_maskL2);
		spz=(i+h_Lsdb2) & h_maskL3;
		smz=(i-h_Lsdb2) & h_maskL3;
		
		e=a[i]*(an[spx]*Jpx[i]+
				an[spy]*Jpy[i]+
				an[spz]*Jpz[i]+
				an[smx]*Jmx[i]+
				an[smy]*Jmy[i]+
				an[smz]*Jmz[i])+
		b[i]*(bn[spx]*Jpx[i]+
			  bn[spy]*Jpy[i]+
			  bn[spz]*Jpz[i]+
			  bn[smx]*Jmx[i]+
			  bn[smy]*Jmy[i]+
			  bn[smz]*Jmz[i])+
		c[i]*(cn[spx]*Jpx[i]+
			  cn[spy]*Jpy[i]+
			  cn[spz]*Jpz[i]+
			  cn[smx]*Jmx[i]+
			  cn[smy]*Jmy[i]+
			  cn[smz]*Jmz[i]);
		energy[0]+=e;
		
	}
	
}


//-------------------------------------------------------------------------------------------------------
// ComputeRandomGPU
//-------------------------------------------------------------------------------------------------------
void ComputeRandomGPU(int numBlocksRand, int numThreadsPerBlockRand, int mode, REAL* array){
	
	seedRand(0);
	
	MY_CUDA_CHECK( cudaThreadSynchronize() );
	
	dim3 dimGrid(numBlocksRand,1,1);
	
	dim3 dimBlock(numThreadsPerBlockRand,1,1);
	
	RandomGPU<<<dimGrid, numThreadsPerBlockRand>>>(array,mode);
	
	//cutilCheckMsg("RandomGPU() execution failed\n");
	
	MY_CUDA_CHECK( cudaThreadSynchronize() );
	
}


//-------------------------------------------------------------------------------------------------------
// FreeRand
//-------------------------------------------------------------------------------------------------------
void FreeRand(void){
	
	MY_CUDA_CHECK( cudaFree(d_aRand) );
	MY_CUDA_CHECK( cudaFree(d_bRand) );
	MY_CUDA_CHECK( cudaFree(d_cRand) );
	MY_CUDA_CHECK( cudaFree(d_Rand) );
	
}


//-------------------------------------------------------------------------------------------------------
// initRand
//-------------------------------------------------------------------------------------------------------
int initRand(int randomNumber){
	
	
	int n_per_rng;
	// int randNum; DA DISCOMMENTARE DOPO PROVA MESSO GLOBAL
	
	n_per_rng = iAlignUp(iDivUp(randomNumber, MT_RNG_COUNT), 2);
	
	randNum = MT_RNG_COUNT * n_per_rng;
	
  //printf("randomNumber: %d n_per_rng: %d MT_RNG_COUNT: %d randNum: %d\n",randomNumber,n_per_rng,MT_RNG_COUNT,randNum);
  
	MY_CUDA_CHECK( cudaMalloc((void **)&d_aRand, randNum * sizeof(REAL)) );
	MY_CUDA_CHECK( cudaMalloc((void **)&d_bRand, randNum * sizeof(REAL)) );
	MY_CUDA_CHECK( cudaMalloc((void **)&d_cRand, randNum * sizeof(REAL)) );
	MY_CUDA_CHECK( cudaMalloc((void **)&d_Rand, randNum * sizeof(REAL)) );
	
	MY_CUDA_CHECK( cudaMemcpyToSymbol("NPerRng",&n_per_rng,sizeof(int),0,cudaMemcpyHostToDevice) );
	
	return randNum;
	
}


//-------------------------------------------------------------------------------------------------------
// seedRand
//-------------------------------------------------------------------------------------------------------
void seedRand(int seed){
	
	if(seed==0) seed=(int)(drand48()*1023); // UNA SCHIFEZZA
	
	static mt_struct_stripped h_MT[MT_RNG_COUNT];
	static mt_struct MT[MT_RNG_COUNT];
	
	FILE *fd = Fopen("./data/MersenneTwister.raw", "rb");
	
	for (int i = 0; i < MT_RNG_COUNT; i++){
		//Inline structure size for compatibility,
		//since pointer types are 8-byte on 64-bit systems (unused *state variable)
		Fread(MT + i, 16 /* sizeof(mt_struct) */ * sizeof(int), 1, fd);
		
	}
	
	fclose(fd);
	
	fd = Fopen("./data/MersenneTwister.dat", "rb");
	
	Fread(h_MT, sizeof(h_MT), 1, fd);
	
	fclose(fd);
	
	//Need to be thread-safe
	mt_struct_stripped *MT_temp = (mt_struct_stripped *)Malloc(MT_RNG_COUNT * sizeof(mt_struct_stripped));
	
	for(int i = 0; i < MT_RNG_COUNT; i++){
		MT_temp[i] = h_MT[i];
		MT_temp[i].seed = seed;
	}
	
	MY_CUDA_CHECK( cudaMemcpyToSymbol("ds_MT",MT_temp,sizeof(h_MT),0,cudaMemcpyHostToDevice) );
	
	free(MT_temp);
	
}


//-------------------------------------------------------------------------------------------------------
// printRand
//-------------------------------------------------------------------------------------------------------
void printRand(int randomNumber,REAL* array,char* namefile){
	
	REAL * h_Rand  = (REAL *)malloc(randomNumber * sizeof(REAL));
	
	MY_CUDA_CHECK( cudaMemcpy(h_Rand, array, randomNumber * sizeof(REAL), cudaMemcpyDeviceToHost) );
	
	FILE* out = Fopen(namefile, "w");
	
	for(int x=0;x<randomNumber;x++) { fprintf(out,"%f\n",h_Rand[x]); }
	
	fclose(out);
	
}


////////////////////////////////////////////////////////////////////////////////
// Write MT_RNG_COUNT vertical lanes of NPerRng random numbers to *d_Random.
// For coalesced global writes MT_RNG_COUNT should be a multiple of warp size.
// Initial states for each generator are the same, since the states are
// initialized from the global seed. In order to improve distribution properties
// on small NPerRng supply dedicated (local) seed to each twister.
// The local seeds, in their turn, can be extracted from global seed
// by means of any simple random number generator, like LCG.
////////////////////////////////////////////////////////////////////////////////
__global__ void RandomGPU( REAL *d_Random, int mode ){
	
	const int tid = blockDim.x * blockIdx.x + threadIdx.x;
	const int THREAD_N = blockDim.x * gridDim.x;
	
	int iState, iState1, iStateM, iOut;
	unsigned int mti, mti1, mtiM, x;
	unsigned int mt[MT_NN];
	
	REAL temp_rand;
	
	for(int iRng = tid; iRng < MT_RNG_COUNT; iRng += THREAD_N){
		//Load bit-vector Mersenne Twister parameters
		mt_struct_stripped config = ds_MT[iRng];
		
		//Initialize current state
		mt[0] = config.seed;
		for(iState = 1; iState < MT_NN; iState++)
		mt[iState] = (1812433253U * (mt[iState - 1] ^ (mt[iState - 1] >> 30)) + iState) & MT_WMASK;
		
		iState = 0;
		mti1 = mt[0];
		
		for(iOut = 0; iOut < NPerRng; iOut++){
			//iState1 = (iState +     1) % MT_NN
			//iStateM = (iState + MT_MM) % MT_NN
			iState1 = iState + 1;
			iStateM = iState + MT_MM;
			if(iState1 >= MT_NN) iState1 -= MT_NN;
			if(iStateM >= MT_NN) iStateM -= MT_NN;
			mti  = mti1;
			mti1 = mt[iState1];
			mtiM = mt[iStateM];
			
			x = (mti & MT_UMASK) | (mti1 & MT_LMASK);
			x = mtiM ^ (x >> 1) ^ ((x & 1) ? config.matrix_a : 0);
			mt[iState] = x;
			iState = iState1;
			
			//Tempering transformation
			x ^= (x >> MT_SHIFT0);
			x ^= (x << MT_SHIFTB) & config.mask_b;
			x ^= (x << MT_SHIFTC) & config.mask_c;
			x ^= (x >> MT_SHIFT1);
			
			temp_rand = ((REAL)x + 1.0f) / 4294967296.0f; //Convert to (0, 1] float
			
			if(mode==ZEROONE) { // write to global memory
				d_Random[iRng + iOut * MT_RNG_COUNT] = temp_rand;
			} else if(mode==ONEONE) {
				d_Random[iRng + iOut * MT_RNG_COUNT] = 2.0f*temp_rand-1.0f;
			}else {
				d_Random[iRng + iOut * MT_RNG_COUNT] = temp_rand - 0.5f;
			}
			
		}
	}
}


//-------------------------------------------------------------------------------------------------------
// HeatBath
//-------------------------------------------------------------------------------------------------------
#if defined(SINGLEHBKERNEL)
__global__ void HeatBath(REAL *a, REAL *b, REAL *c,
                         REAL *an, REAL *bn, REAL *cn,
                         REAL *Jpx, REAL *Jpy, REAL *Jpz,
                         REAL *Jmx, REAL *Jmy, REAL *Jmz,
                         REAL *atry, REAL *btry, REAL *ctry,
                         REAL *random, int even) {
	
	unsigned int iy, iz;
	
	int diff;
	REAL as, bs, cs;
	REAL anew, bnew, cnew;
	REAL hk, bhk, projection, pl, invnorm;
	
	REAL invexpbhk2;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart;
	
	istart=blockDim.x*blockIdx.x;
	
	
	for(i=istart+tid; i<d_nspin; i+=num_threads) {
		
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		//do the same thing if both $iz and $iy are even or odd
		diff=(((iz&1)^(iy&1))^even);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2) & d_maskL3;
		smz=(i-d_Lsdb2) & d_maskL3;
		
		as=an[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
		an[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
		an[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
		an[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
		an[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
		an[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];
		
		bs=bn[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
		bn[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
		bn[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
		bn[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
		bn[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
		bn[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];
		
		cs=cn[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
		cn[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
		cn[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
		cn[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
		cn[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
		cn[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);
		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;
		
	}
	
}

#else


//-------------------------------------------------------------------------------------------------------
// HeatBathRedBnd
//-------------------------------------------------------------------------------------------------------
__global__ void HeatBathRedBnd(REAL *a, REAL *b, REAL *c,
							   REAL *an, REAL *bn, REAL *cn,
							   REAL *Jpx, REAL *Jpy, REAL *Jpz,
							   REAL *Jmx, REAL *Jmy, REAL *Jmz,
							   REAL *atry, REAL *btry, REAL *ctry,
							   REAL *random) {
	
	unsigned int iy, iz;
	
	int diff;
	REAL as, bs, cs;
	REAL anew, bnew, cnew;
	REAL hk, bhk, projection, pl, invnorm;
	
	REAL invexpbhk2;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart;
	
	istart=blockDim.x*blockIdx.x;
	
	
	for(i=istart+tid; i<d_Lsdb2; i+=num_threads) {
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2);// & d_maskL3;
		smz=(i-d_Lsdb2);// & d_maskL3;
		
		as=an[spx]*Jpx[i]+
		an[spy]*Jpy[i]+
		an[spz]*Jpz[i]+
		an[smx]*Jmx[i]+
		an[smy]*Jmy[i]+
		an[smz]*Jmz[i];
		
		bs=bn[spx]*Jpx[i]+
		bn[spy]*Jpy[i]+
		bn[spz]*Jpz[i]+
		bn[smx]*Jmx[i]+
		bn[smy]*Jmy[i]+
		bn[smz]*Jmz[i];
		
		cs=cn[spx]*Jpx[i]+
		cn[spy]*Jpy[i]+
		cn[spz]*Jpz[i]+
		cn[smx]*Jmx[i]+
		cn[smy]*Jmy[i]+
		cn[smz]*Jmz[i];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);
		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;
		
	}
	
	istart=(d_nspin-d_Lsdb2)+blockDim.x*blockIdx.x;
	
	
    for(i=istart+tid; i<d_nspin; i+=num_threads) {
		
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2);// & d_maskL3;
		smz=(i-d_Lsdb2);// & d_maskL3;
		
		as=an[spx]*Jpx[i]+
		an[spy]*Jpy[i]+
		an[spz]*Jpz[i]+
		an[smx]*Jmx[i]+
		an[smy]*Jmy[i]+
		an[smz]*Jmz[i];
		
		bs=bn[spx]*Jpx[i]+
		bn[spy]*Jpy[i]+
		bn[spz]*Jpz[i]+
		bn[smx]*Jmx[i]+
		bn[smy]*Jmy[i]+
		bn[smz]*Jmz[i];
		
		cs=cn[spx]*Jpx[i]+
		cn[spy]*Jpy[i]+
		cn[spz]*Jpz[i]+
		cn[smx]*Jmx[i]+
		cn[smy]*Jmy[i]+
		cn[smz]*Jmz[i];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);
		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;
		
	}
	
}


//-------------------------------------------------------------------------------------------------------
// HeatBathRedBulk
//-------------------------------------------------------------------------------------------------------
__global__ void HeatBathRedBulk(REAL *a, REAL *b, REAL *c,
								REAL *an, REAL *bn, REAL *cn,
								REAL *Jpx, REAL *Jpy, REAL *Jpz,
								REAL *Jmx, REAL *Jmy, REAL *Jmz,
								REAL *atry, REAL *btry, REAL *ctry,
								REAL *random) {
	
	unsigned int iy, iz;
	
	int diff;
	REAL as, bs, cs;
	REAL anew, bnew, cnew;
	REAL hk, bhk, projection, pl, invnorm;
	
	REAL invexpbhk2;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart, iend;
	
	istart=blockDim.x*blockIdx.x+d_Lsdb2; //start on second plane (by using d_Lsdb2)
										  // since first plane is for the boundary 
	iend=d_nspin-d_Lsdb2;
	
	for(i=istart+tid; i<iend; i+=num_threads) {
    
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];

		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2);// & d_maskL3;
		smz=(i-d_Lsdb2);// & d_maskL3;
		
		as=an[spx]*Jpx[i]+
		an[spy]*Jpy[i]+
		an[spz]*Jpz[i]+
		an[smx]*Jmx[i]+
		an[smy]*Jmy[i]+
		an[smz]*Jmz[i];
		
		bs=bn[spx]*Jpx[i]+
		bn[spy]*Jpy[i]+
		bn[spz]*Jpz[i]+
		bn[smx]*Jmx[i]+
		bn[smy]*Jmy[i]+
		bn[smz]*Jmz[i];
		
		cs=cn[spx]*Jpx[i]+
		cn[spy]*Jpy[i]+
		cn[spz]*Jpz[i]+
		cn[smx]*Jmx[i]+
		cn[smy]*Jmy[i]+
		cn[smz]*Jmz[i];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);

		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;

		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;
		
		
	}

}



//-------------------------------------------------------------------------------------------------------
// HeatBathBlueBnd
//-------------------------------------------------------------------------------------------------------
__global__ void HeatBathBlueBnd(REAL *a, REAL *b, REAL *c,
								REAL *an, REAL *bn, REAL *cn,
								REAL *Jpx, REAL *Jpy, REAL *Jpz,
								REAL *Jmx, REAL *Jmy, REAL *Jmz,
								REAL *atry, REAL *btry, REAL *ctry,
								REAL *random) {
	
	unsigned int iy, iz;
	
	int diff;
	REAL as, bs, cs;
	REAL anew, bnew, cnew;
	REAL hk, bhk, projection, pl, invnorm;
	
	REAL invexpbhk2;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart;
	
	istart=blockDim.x*blockIdx.x;
	
	
	for(i=istart+tid; i<d_Lsdb2; i+=num_threads) {
		
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2);// & d_maskL3;
		smz=(i-d_Lsdb2);// & d_maskL3;
		
		as=an[spx]*Jpx[spx]+
		an[spy]*Jpy[spy]+
		an[spz]*Jpz[spz]+
		an[smx]*Jmx[smx]+
		an[smy]*Jmy[smy]+
		an[smz]*Jmz[smz];
		
		bs=bn[spx]*Jpx[spx]+
		bn[spy]*Jpy[spy]+
		bn[spz]*Jpz[spz]+
		bn[smx]*Jmx[smx]+
		bn[smy]*Jmy[smy]+
		bn[smz]*Jmz[smz];
		
		
		cs=cn[spx]*Jpx[spx]+
		cn[spy]*Jpy[spy]+
		cn[spz]*Jpz[spz]+
		cn[smx]*Jmx[smx]+
		cn[smy]*Jmy[smy]+
		cn[smz]*Jmz[smz];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);
		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;
		
	}
	
	istart=(d_nspin-d_Lsdb2)+blockDim.x*blockIdx.x;
	
	
	for(i=istart+tid; i<d_nspin; i+=num_threads) {
		
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2);// & d_maskL3;
		smz=(i-d_Lsdb2);// & d_maskL3;
		
		as=an[spx]*Jpx[spx]+
		an[spy]*Jpy[spy]+
		an[spz]*Jpz[spz]+
		an[smx]*Jmx[smx]+
		an[smy]*Jmy[smy]+
		an[smz]*Jmz[smz];
		
		bs=bn[spx]*Jpx[spx]+
		bn[spy]*Jpy[spy]+
		bn[spz]*Jpz[spz]+
		bn[smx]*Jmx[smx]+
		bn[smy]*Jmy[smy]+
		bn[smz]*Jmz[smz];
		
		
		cs=cn[spx]*Jpx[spx]+
		cn[spy]*Jpy[spy]+
		cn[spz]*Jpz[spz]+
		cn[smx]*Jmx[smx]+
		cn[smy]*Jmy[smy]+
		cn[smz]*Jmz[smz];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);
		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;
		
	}
	
	
}



//-------------------------------------------------------------------------------------------------------
// HeatBathBlueBulk
//-------------------------------------------------------------------------------------------------------
__global__ void HeatBathBlueBulk(REAL *a, REAL *b, REAL *c,
								 REAL *an, REAL *bn, REAL *cn,
								 REAL *Jpx, REAL *Jpy, REAL *Jpz,
								 REAL *Jmx, REAL *Jmy, REAL *Jmz,
								 REAL *atry, REAL *btry, REAL *ctry,
								 REAL *random) {
	
	
	unsigned int iy, iz;
	
	int diff;
	REAL as, bs, cs;
	REAL anew, bnew, cnew;
	REAL hk, bhk, projection, pl, invnorm;
	
	REAL invexpbhk2;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart, iend;
	
	istart=blockDim.x*blockIdx.x;
	//istart=blockDim.x*blockIdx.x+d_Lsdb2; //start on second plane (by using d_Lsdb2)
										  // since first plane is for the boundary 

	iend=d_nspin-d_Lsdb2;
	
	for(i=istart+tid; i<iend; i+=num_threads) {
		
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2);// & d_maskL3;
		smz=(i-d_Lsdb2);// & d_maskL3;
		
		as=an[spx]*Jpx[spx]+
		an[spy]*Jpy[spy]+
		an[spz]*Jpz[spz]+
		an[smx]*Jmx[smx]+
		an[smy]*Jmy[smy]+
		an[smz]*Jmz[smz];
		
		bs=bn[spx]*Jpx[spx]+
		bn[spy]*Jpy[spy]+
		bn[spz]*Jpz[spz]+
		bn[smx]*Jmx[smx]+
		bn[smy]*Jmy[smy]+
		bn[smz]*Jmz[smz];
		
		
		cs=cn[spx]*Jpx[spx]+
		cn[spy]*Jpy[spy]+
		cn[spz]*Jpz[spz]+
		cn[smx]*Jmx[smx]+
		cn[smy]*Jmy[smy]+
		cn[smz]*Jmz[smz];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);
		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		/*a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;*/
		
	}
	
}


#endif

struct mpiexch *SetUpMpi(int id, int np, int side) {
	struct mpiexch *p;
	p=(struct mpiexch*)Malloc(sizeof(struct mpiexch));
	MY_CUDA_CHECK( cudaMallocHost((void **)&p->upsend,side*side*sizeof(MYREAL)*NCOMP) );
	MY_CUDA_CHECK( cudaMallocHost((void **)&p->uprecv,side*side*sizeof(MYREAL)*NCOMP) );
	MY_CUDA_CHECK( cudaMallocHost((void **)&p->downsend,side*side*sizeof(MYREAL)*NCOMP) );
	MY_CUDA_CHECK( cudaMallocHost((void **)&p->downrecv,side*side*sizeof(MYREAL)*NCOMP) );
	p->up=id+1<np?id+1:0;
	p->down=id-1<0?np-1:id-1;
	p->nbyte=side*side*sizeof(MYREAL)*NCOMP;
	p->post=FALSE;
	return p;
}

